1_PLSC_grad_behav
Ju-Chi.Yu
3/25/2022
Grabbing SPINS gradients
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Warning: package 'ggseg' was built under R version 4.1.3## Loading required package: prettyGraphs## Loading required package: ExPositionRead in the SPINS big table
## New names:
## * `` -> ...1## Rows: 337120 Columns: 16## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (5): ROI, Task, Subject, Network, Site
## dbl (11): ...1, grad1, grad2, grad3, grad4, grad5, grad6, grad7, grad8, grad...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.## New names:
## * `` -> ...1## Rows: 337120 Columns: 16## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (5): ROI, Task, Subject, Network, Site
## dbl (11): ...1, grad1, grad2, grad3, grad4, grad5, grad6, grad7, grad8, grad...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.read subject data
## Rows: 354 Columns: 38## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (7): record_id, site, scanner, diagnostic_group, demo_sex, subject, ses...
## dbl (22): demo_age_study_entry, scog_tasit1_total, scog_tasit2_total, scog_t...
## lgl (9): term_early_withdraw, has_RS_grads, has_EA_grads, RS_QC_exclude, EA...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.## [1] "record_id" "site"
## [3] "scanner" "diagnostic_group"
## [5] "demo_sex" "demo_age_study_entry"
## [7] "term_early_withdraw" "scog_tasit1_total"
## [9] "scog_tasit2_total" "scog_tasit3_total"
## [11] "np_composite_tscore" "np_domain_tscore_process_speed"
## [13] "np_domain_tscore_att_vigilance" "np_domain_tscore_work_mem"
## [15] "np_domain_tscore_verbal_learning" "np_domain_tscore_visual_learning"
## [17] "np_domain_tscore_reasoning_ps" "np_domain_tscore_social_cog"
## [19] "bprs_factor_pos_symp" "bprs_factor_total"
## [21] "bprs_factor_neg_symp" "qls_total"
## [23] "sans_total_sc" "bsfs_total"
## [25] "subject" "session"
## [27] "n_scans_task-emp" "n_scans_task-rest"
## [29] "fd_mean_task-emp" "fd_mean_task-rest"
## [31] "has_RS_grads" "has_EA_grads"
## [33] "RS_QC_exclude" "EA_QC_exclude"
## [35] "RS_mo_exclude" "EA_mo_exclude"
## [37] "withdrawn_exclude" "exclude_MRI"## [1] "scog_tasit1_total" "scog_tasit2_total"
## [3] "scog_tasit3_total" "np_composite_tscore"
## [5] "np_domain_tscore_process_speed" "np_domain_tscore_att_vigilance"
## [7] "np_domain_tscore_work_mem" "np_domain_tscore_verbal_learning"
## [9] "np_domain_tscore_visual_learning" "np_domain_tscore_reasoning_ps"
## [11] "np_domain_tscore_social_cog" "bprs_factor_pos_symp"
## [13] "bprs_factor_total" "bprs_factor_neg_symp"
## [15] "qls_total" "sans_total_sc"
## [17] "bsfs_total"Check subject overlap
grad.sub <- spins_grads_wide$Subject[order(spins_grads_wide$Subject)]
behav.sub <- lol_spins_behav$subject[order(lol_spins_behav$subject)]
behav.sub[behav.sub %in% grad.sub == FALSE]## [1] "sub-CMP0180" "sub-CMP0182" "sub-CMP0196" "sub-CMP0198" "sub-CMP0213"
## [6] "sub-ZHH0034"grad.sub[grad.sub %in% behav.sub == FALSE]## [1] "sub-CMH0009" "sub-CMH0036" "sub-CMH0046" "sub-CMH0047" "sub-CMH0048"
## [6] "sub-CMH0050" "sub-CMH0052" "sub-CMH0059" "sub-CMH0070" "sub-CMH0072"
## [11] "sub-CMH0097" "sub-CMH0104" "sub-CMH0114" "sub-CMH0118" "sub-CMH0126"
## [16] "sub-CMH0136" "sub-CMH0155" "sub-CMH0176" "sub-CMH0197" "sub-CMP0175"
## [21] "sub-CMP0178" "sub-CMP0183" "sub-CMP0202" "sub-MRC0001" "sub-MRC0008"
## [26] "sub-MRC0010" "sub-MRC0012" "sub-MRC0013" "sub-MRC0023" "sub-MRC0042"
## [31] "sub-MRC0044" "sub-MRC0046" "sub-MRC0047" "sub-MRC0051" "sub-MRC0058"
## [36] "sub-MRC0060" "sub-MRC0070" "sub-MRC0071" "sub-MRP0019" "sub-MRP0067"
## [41] "sub-MRP0079" "sub-MRP0093" "sub-MRP0099" "sub-MRP0105" "sub-MRP0116"
## [46] "sub-MRP0121" "sub-MRP0132" "sub-MRP0143" "sub-MRP0165" "sub-ZHH0014"
## [51] "sub-ZHH0017" "sub-ZHH0019" "sub-ZHH0020" "sub-ZHH0022" "sub-ZHH0023"
## [56] "sub-ZHH0033" "sub-ZHH0037" "sub-ZHH0040" "sub-ZHH0044" "sub-ZHH0045"
## [61] "sub-ZHH0047" "sub-ZHH0048" "sub-ZHH0052" "sub-ZHH0053" "sub-ZHH0054"
## [66] "sub-ZHH0057" "sub-ZHP0063" "sub-ZHP0065" "sub-ZHP0085" "sub-ZHP0093"
## [71] "sub-ZHP0094" "sub-ZHP0095" "sub-ZHP0097" "sub-ZHP0098" "sub-ZHP0100"
## [76] "sub-ZHP0103" "sub-ZHP0105" "sub-ZHP0119" "sub-ZHP0129" "sub-ZHP0133"
## [81] "sub-ZHP0134" "sub-ZHP0141" "sub-ZHP0142" "sub-ZHP0143" "sub-ZHP0144"
## [86] "sub-ZHP0155" "sub-ZHP0162" "sub-ZHP0168"## kept subjects
kept.sub <- behav.sub[behav.sub %in% grad.sub] # 342
kept.sub[complete.cases(spins_behav_num[kept.sub,c(1:11, 17)]) == FALSE]## [1] "sub-CMH0013" "sub-ZHP0081" "sub-ZHP0083" "sub-ZHP0118"kept.sub <- kept.sub[complete.cases(spins_behav_num[kept.sub,c(1:11, 17)])] # 338
## grab the matching data
behav.dat <- spins_behav_num[kept.sub,c(1:11, 17)]
spins_grads_wide_org <- spins_grads_wide[,-1]
rownames(spins_grads_wide_org) <- spins_grads_wide$Subject
grad.dat <- spins_grads_wide_org[kept.sub,]
rm(spins_grads_wide, spins_grads_wide_org, lol_spins_behav, spins_behav_num, spins_grads_num)grab some network colours
networks <- read_delim("../networks.txt",
"\t", escape_double = FALSE, trim_ws = TRUE) %>%
select(NETWORK, NETWORKKEY, RED, GREEN, BLUE, ALPHA) %>%
distinct() %>%
add_row(NETWORK = "Subcortical", NETWORKKEY = 13, RED = 0, GREEN=0, BLUE=0, ALPHA=255) %>%
mutate(hex = rgb(RED, GREEN, BLUE, maxColorValue = 255)) %>%
arrange(NETWORKKEY)## Rows: 718 Columns: 12## -- Column specification --------------------------------------------------------
## Delimiter: "\t"
## chr (4): LABEL, HEMISPHERE, NETWORK, GLASSERLABELNAME
## dbl (8): INDEX, KEYVALUE, RED, GREEN, BLUE, ALPHA, NETWORKKEY, NETWORKSORTED...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.networks## # A tibble: 13 x 7
## NETWORK NETWORKKEY RED GREEN BLUE ALPHA hex
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Visual1 1 0 0 255 255 #0000FF
## 2 Visual2 2 100 0 255 255 #6400FF
## 3 Somatomotor 3 0 255 255 255 #00FFFF
## 4 Cingulo-Opercular 4 153 0 153 255 #990099
## 5 Dorsal-Attention 5 0 255 0 255 #00FF00
## 6 Language 6 0 155 155 255 #009B9B
## 7 Frontoparietal 7 255 255 0 255 #FFFF00
## 8 Auditory 8 250 62 251 255 #FA3EFB
## 9 Default 9 255 0 0 255 #FF0000
## 10 Posterior-Multimodal 10 177 89 40 255 #B15928
## 11 Ventral-Multimodal 11 255 157 0 255 #FF9D00
## 12 Orbito-Affective 12 65 125 0 168 #417D00
## 13 Subcortical 13 0 0 0 255 #000000get row and column designs
## match ROIs to networks
ROI.network.match <- cbind(spins_grads$ROI, spins_grads$Network) %>% unique
ROI.idx <- ROI.network.match[,2]
names(ROI.idx) <- ROI.network.match[,1]
### match networks with colors
net.col.idx <- networks$hex
names(net.col.idx) <- networks$NETWORK
## design matrix for subjects
sub.dx <- spins_dx_org[kept.sub,]
diagnostic.col <- sub.dx$diagnostic_group %>% as.matrix %>% makeNominalData() %>% createColorVectorsByDesign()
rownames(diagnostic.col$gc) <- sub(".","", rownames(diagnostic.col$gc))
## design matrix for columns - behavioral
behav.dx <- matrix(nrow = ncol(behav.dat), ncol = 1, dimnames = list(colnames(behav.dat), "type")) %>% as.data.frame
behav.col <- c("scog" = "#F28E2B",
"np" = "#59A14F",
"bsfs" = "#D37295")
behav.dx$type <- sub("(^[^_]+).*", "\\1", colnames(behav.dat))
behav.dx$type.col <- recode(behav.dx$type, !!!behav.col)
## design matrix for columns - gradient
grad.dx <- matrix(nrow = ncol(grad.dat), ncol = 4, dimnames = list(colnames(grad.dat), c("gradient", "ROI", "network", "network.col"))) %>% as.data.frame
grad.dx$gradient <- sub("(^[^.]+).*", "\\1", colnames(grad.dat))
grad.dx$ROI <- sub("^[^.]+.(*)", "\\1", colnames(grad.dat))
grad.dx$network <- recode(grad.dx$ROI, !!!ROI.idx)
grad.dx$network.col <- recode(grad.dx$network, !!!net.col.idx)
## get different alpha for gradients
grad.col.idx <- c("grad1" = "grey30",
"grad2" = "grey60",
"grad3" = "grey90")
grad.dx$gradient.col <- recode(grad.dx$gradient, !!!grad.col.idx)Run PLS-C
pls.res <- tepPLS(behav.dat, grad.dat, DESIGN = sub.dx$diagnostic_group, make_design_nominal = TRUE, graphs = FALSE)
## swith direction for dimension 3
pls.res$TExPosition.Data$fi[,3] <- pls.res$TExPosition.Data$fi[,3]*-1
pls.res$TExPosition.Data$fj[,3] <- pls.res$TExPosition.Data$fj[,3]*-1
pls.res$TExPosition.Data$pdq$p[,3] <- pls.res$TExPosition.Data$pdq$p[,3]*-1
pls.res$TExPosition.Data$pdq$q[,3] <- pls.res$TExPosition.Data$pdq$q[,3]*-1
pls.res$TExPosition.Data$lx[,3] <- pls.res$TExPosition.Data$lx[,3]*-1
pls.res$TExPosition.Data$ly[,3] <- pls.res$TExPosition.Data$ly[,3]*-1
PlotScree(pls.res$TExPosition.Data$eigs)
Results
Dimension 1
lxly.out[[1]]
gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)
PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,1],
threshold = 0,
color4bar = behav.dx$type.col,
horizontal = FALSE, main = "Loadings - behavioural")## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
Dimension 2
lxly.out[[2]]
gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)
PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,2],
threshold = 0, color4bar = behav.dx$type.col, horizontal = FALSE, main = "Loadings - behavioural")
Dimension 3
lxly.out[[3]]
gridExtra::grid.arrange(bar.grad1, bar.grad2, bar.grad3, ncol = 1)
PrettyBarPlot2(pls.res$TExPosition.Data$pdq$p[,3],
threshold = 0, color4bar = behav.dx$type.col, horizontal = FALSE, main = "Loadings - behavioural")
Back into the brain
Dimension 1
## merging atlas and data by 'label'
## merging atlas and data by 'label'
Dimension 2
## merging atlas and data by 'label'
## merging atlas and data by 'label'
Dimension 3
## merging atlas and data by 'label'
## merging atlas and data by 'label'
Fancy figures
3D plot of the gradients
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutdata4plotly <- data.frame(spins_mean, grad4grp[rownames(spins_mean),])
color4plotly <- net.col.idx
fig <- plot_ly(data4plotly, x = ~grad1, y = ~grad2, z = ~grad3, color = ~Network, colors = color4plotly, alpha = 0.5, size = 5)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = "Gradient 1"),
yaxis = list(title = "Gradient 2"),
zaxis = list(title = "Gradient 3")))
figDimension 1
Dimension 2
Dimension 3
